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Abstract 

To efficiently implement many-qubit gates for use in quantum simulations on 
quantum computers we develop and present methods reexpressing exp[— i[H% + 
H2 + . . -)At] as a product of factors exp[— iHiAt], exp[— iH^At], . . . which is ac- 
curate to 3rd or 4th order in At. The methods we derive are an extended form of 
symplectic method and can also be used for the integration of classical Hamiltoni- 
ans on classical computers. We derive both integral and irrational methods, and 
find the most efficient methods in both cases. 



1 Introduction 



Quantum computers have generated much interest recently, largely due to the result by 
Shor [Q] that they can factor integers in polynomial time. 

In a quantum computer the analog of a logical bit is the qubit. The canonical example 
of a qubit is a quantum spin. A quantum spin consists of two states, so a set of n spins 
gives the quantum computer a 2 n -dimensional Hilbert space. 

To perform a calculation, one initializes the qubits, and then applies unitary logical 
gates to the qubits. Unitary logical gates are realised in different fashions depending 
on the quantum computer hardware, but they are all represented mathematically by a 
Hamiltonian acting on the quantum state of the qubits. In a typical quantum computer, 
technology restricts the Hamiltonian to act on a small number of qubits at a time, maybe 
two or three. A calculation is then built up of two- or three-qubit Hamiltonians, or gates, 
acting sequentially on the qubits. 

An important and difficult to realise requirement is that the qubits maintain their 
coherence throughout an entire calculation. Maintaining coherence in quantum com- 
puters is a problem which has led to the development of error correcting codes (see || 
and included references). These codes are possible due to the fact that one does not 
need to know the state of a qubit in order to tell whether an error has occurred. With 
some ingenuity, it is possible to determine what kinds of errors have occurred during the 
course of a calculation and to correct the errors as the calculation proceeds. Simple error 
correction codes have already been shown to work on small numbers of qubits ||. 

Effort has also been put into developing algorithms which make use of the quantum 
computer's power. Shor's algorithm showed that quantum computers are more powerful 
than classical computers, since integers cannot be factored in polynomial time on a 
classical computer, whereas they can on a quantum computer. Grover has also devised 
a method for searching a database in time proportional to the square root of the number 
of items involved in the search Q. 

In addition to research into effective algorithms for use on quantum computers, sim- 
ulations of quantum systems have also been shown to be possible in polynomial time 0. 
Indeed, this was the first area for which it was proposed that quantum computers could 
fundamentally be more powerful (i.e. much faster) than classical computers 

This paper focuses on a problem which concerns simulational issues in quantum 
computation. Essentially, we have developed methods for reexpressing exp[— i(Hi + 
H 2 + . . .) At] as a product of factors exp[— iHiAt], exp[— iH 2 At], . . . which is accurate to 
3rd or 4th order in At, as mentioned in the abstract. 

A simulation on a quantum computer consists of applying an operator exp(-iHt) on 
a set of qubits, where H, the Hamiltonian of the system of interest is suitably encoded 
(and discretized) to act on the set of qubits. For many body systems, if is a sum of 
terms. For instance, in a one-dimensional Ising spin model, the Hamiltonian is 

N 

H = °n ■ On+1 (1) 

71=1 
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where N is the number of spins. Another example is the Hubbard model Hamiltonian, 
used in the study of high-T c superconductivity, which can be written [|J as the sum 



m 



(2) 



where Vq is the strength of the potential, and n ia is the operator for the number of 
fermions of spin a at site i. In the second (kinetic energy) term, the sum indicates all 
neighboring pairs of sites, t is the strength of the "hopping", and c ia , c* a are annihilation 
and creation operators, respectively, of a fermion at site i and spin a. 

These models give examples in which a large simulation on a classical computer is 
impossible due to the exponential increase in the size of the Hilbert space of the quantum 
system with the number of lattice sites. A many-particle system can sometimes be 
simulated with fewer qubits in first-quantized form O, but in either case the Hamiltonian 
if is a sum of terms, so our methods are equally applicable to both cases. 

If the quantum computer cannot act on all spins at once, as is the case for quantum 
gate arrays || , it becomes necessary to find ways of approximating the application of the 
above Hamiltonians with few-qubit gates. To second order, for instance, we find that 

( -iHiAt -iH 2 At -iH N At\ ( -iH N At -iU 2 At -itfiAA __ -i2(H 1 +H 2 +...H N )At+o[(Atf] 



where H n are two-qubit gates (e.g. a n ■ cr n+ i). 

Below, we analyze the problem of deriving higher order methods of this type, and 
find a set of equations which, once solved, give 3rd and 4th order methods analogous 
to the above second-order method. We solve and present the formulae for 3rd and 4th 
order methods as well as developing methods for approximating expressions involving 
commutators, exp([v4, B}). Since there is a large set of solutions to our equations, we 
spend some effort trying to isolate and present only the most efficient methods. 

After presenting our methods, we then provide results from a simple application to 
give the reader confidence that our methods are correct. 

This kind of method has been investigated elsewhere, for different reasons, in the 
context of Hamiltonian systems under the name 'symplectic' method. In the section on 
symplectic methods, we comment on what we have done differently from other inves- 
tigations of symplectic methods, and why our methods are applicable to more general 
problems. We then present a summary of our results in the conclusions section. 

We also provide appendices with useful expressions used in the derivation of our 
results, and some proofs of statements in the text. 




(3) 
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2 Mathematical Analysis and Equations 

We want to express exp fe^Li Ai) as a product of individual exp (A n ) 's. In order to 
do this, we use the Campbell-Baker-Hausdorff formula. The Campbell-Baker-Hausdorff 
formula to 5th order is 

1 o . 1 , , . . , 1 



a (Ai + A 2 ) + ^-a 2 A 12 + —a 3 (A 112 + A 221 ) + —a 4 A 



— — ±Zi , — ^ — ±x.a 1 — ^^j./ 1 2^ — — 1221 



exp (aAi) exp (aA 2 ) = exp 

—a 5 (Aim — 2A 21112 — 6A1221 — QA 22 u 2 — 2A 1222 i + A 2222 i) + (9 (a 6 ^) 



(4) 



720 
where 

Akl...mn = [A; [A; • • • [An; A] • • •]] (5) 

To find combinations of operators exp A which approximate exp fe^Li A J to some 
order it is first necessary to choose a strategy for searching among the large number 
of possible combinations. First of all, we cannot search brute force since there are too 
many possible combinations, and, in any case, this would not give us a formula valid 
for all N. Therefore, we pick a fundamental ordering of the product of exponentials 
with parameters allowing for transposes of the entire product as well as raising all the 
exponentials in the fundamental unit to the same power. 

By iterating the Campbell-Baker-Hausdorff formula, we can get an expression for the 
fundamental unit in terms of a single exponential 

oo 

(e aM e aM . . . e aAN ) a = exp £ aa p B p N (6) 

P =i 

which defines the B P N in terms of the A n . Here, p is an exponent on a, and a label on 
the matrices B P N . We take a = ±1. 

Now combine a succession i = 1, . . . , / of fundamental units with parameters dj and 
ad- Again iterating Campbell-Baker-Hausdorff gives 

exp (^T aialB p N j ... exp (^T a^B^j = exp ^ of B^j (7) 

The B* are generated from the B P N by commutation. X represents a label pq . . .rs where 

B p N ^ s ^[B p N ,{B%,...[B^B^ N ]...]] (8) 

B p ^'" rs is of order p + q + . . . +r + s. Up to 5th order we can take 

X E {1; 2; 3, 12; 4, 13, 112; 5, 14, 23, 113, 221, 1112} (9) 

These B§ span the space of commutators of the B p N 's to 5th order and for N > 2 they 
are independent. Formulae for the B* in terms of A and A 2 are given in Appendix A. 2. 
The of are defined in terms of a, and by Eq. (fl). Here again, the A's are labels. 
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After some calculation, the Campbell-Baker-Hausdorff formula, Eq. (^), then gives 



J2 a i a i P 



(10) 



for p — 1, 



a 



pq 



4*4 + 5& 



cr- 



(11) 



for pq = 12, 13, 14, 23, 



ppq 



-\a\af - \ (a*) 



(cr 



P\3 



(12) 



for ppq = 112,113,221 Q, 



^1112 
I 



,1„H2 



-ova- 



rii 



^;)V-^!)^+^x>,[(*,' 



(13) 



For approximations to exp fe^Li AiJ , w e require all = except for a} which is 
the coefficient of B]^ = Y,n=i and which should be greater than zero. 

An interesting feature of 3rd order methods is that they require inverses, i.e. they 
require backward time evolution during part of the method.0 This can be proved using 
Eq. ( |T0"D with p = 3. It has no nontrivial solutions when the product ajOj is positive for 
all i. Therefore for 3rd order methods, a^aj must be negative for at least one i. From the 
left hand side of Eq. (|]), we see that this means that there must be at least one inverse. 
Similarly, from Eqs. (|ToD with p = 3 and p — 4 it can be proved that 4th order methods 
must have at least two inverses. 

In Appendix A.l, we also prove the fact that for integral solutions o\ must be a 
multiple of 2 for a 2nd order method, a multiple of 6 for a 3rd or 4th order method, and 
a multiple of 30 for a 5th order method. Our searches suggest that the constraints on a} 
may actually be stronger; all 4th order methods that we have found have a} a multiple 
of 12, and we have not been able to find any 5th order methods. 

In Section [7|, we will consider approximation to expfAi,^] for which we require all 
of = except for erf. 



1 For the purposes of calculating erf 21 , note that erf 1 = — a] 2 . 

2 After this work was completed, we became aware that this point had also been noted in 
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3 Numerical Method for Solution of the Equations 



We solve our equations for both integer and irrational solutions, using different methods 
for each search. 

Our method to solve Eqs. ([10| - |T2|) for integers is to pick values of and cij and see 
if they satisfy the equations. To do this we restrict the number of fundamental units by 
fixing I. We also restrict the range of the a^s. 

We start with Eq. fllPf) , since, in this equation, order with respect to i does not 
matter. So, for a given set of values, we need to consider only one permutation, not all 
permutations of the values. This greatly reduces the size of the search. 

Furthermore, we start by considering p = 1 and 3, since it is only the sign of a^ai that 
matters in these equations. This means we can consider only the sign of the combination 
ajdj, and not the signs of «j and individually. This reduces the search further. These 
equations are particularly restrictive for the case of few inverses. 

After solving the p = 1 and 3 equations, we introduce separate signs for the a^s and 
aj's and solve the equation with p = 2, and p = 4 for the 4th order case. 

Finally, into the restricted set of solutions to Eq. (|T0| ) we introduce permutations of 
the aj's and a^'s with respect to the index i and solve Eqs. (ITO-O). 



We solved Eqs (|10|) and analytically to find the unique shortest irrational 3rd 
order method. To find 4th order irrational methods, we made a symmetric ansatz and 
solved Eqs. ( Jl0| - |l2|) analytically to find the shortest symmetric irrational 4th order meth- 



ods. We checked numerically, using the globally convergent technique prescribed in ||10|| , 
that these are all the shortest irrational 4th order methods. 
The methods are presented in Section (pj). 



4 Criteria for Selecting Among the Solutions 



With our strategy for finding solutions to Eqs. ( |I~0| - |T2] ) we find a larger number of solutions 
than we can easily present. We need to select solutions to present and we also want to 
present solutions which are in some sense optimal. To do this, we consider the form of 
the operator resulting from a given method 



n( 



(3 iajAi At ^—iajAi At ^—i-o-j^N At \ " J 



exp 



A? 



-iajYl A n At + r(-i At) 0+1 



n=l 



(14) 



where At <C 1 is a time step, o is the order of the method, and 



(15) 



x 



where X e {4, 13, 112} for a 3rd order method and X G {5, 14, 23, 113, 221, 1112} for a 
4th order method. 
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r is an error which takes values in the vector space of the commutators for which we 
do not have a metric. Therefore, we make the ad hoc choice of basis that is given in 
Appendix A. 3. This allows us to replace r by a single real scalar R as is also described 
in Appendix A. 3. The error from the method can then be taken to be 

E = nRAt° +1 (16) 

where n is the number of times we apply the approximate method. 
If the physical time we want to simulate is T p , then 

T p = nDAt (17) 

where D = a} is given by the method. 

The computer time it takes for a given simulation can be written 

T c = nINtg + nLNt s (18) 

where I is the number of fundamental units in the method and N is the number of terms 
in a unit, t g is the time it takes to make the gate change, 

I = E|aj| (19) 

i=l 

so that LN is the total time the gates are applied for in the method, and t s is the time 
each individual gate is applied for. The time an individual gate is applied for will be 
t s = bAt, where b is a proportionality constant dictated by the actual couplings in the 
quantum computer hardware. 

Using Eqs. fll6"D and (|T7|), the computer time can be rewritten 




There are two possible limits to this equation. If At can be made very small (from 
the hardware point of view), then making the error small forces the computer time to 
be dominated by gate switching. In this case, we want the factor 

Z = (I/D)(R/D) l l° (21) 

to be small. 

If there is a lower limit to At = e, and it is reached before the computer time is gate 
switching dominated, then the computer time may be dominated by gate application. In 
this second limit, we want L/D small, and to minimize the error E, we want (R/ D)(At)° 
small. In this limit, each gate can only be applied for an integral number of the minimum 
timestep e. Thus, to use an irrational method, one must approximate the method by 
an integral method containing large integers, and so with a large D. Because At = e is 
fixed, the error E goes like R/D oc D°, and thus is large for irrational methods. We thus 
do not consider irrational methods in this limit. 

To summarize, we want methods with small L/D and R/D. 
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5 3rd and 4th Order Formulae for exp (E^ =1 A n J 

From this analysis, we want to choose methods for which Z, or L/D and R/D are small. 
Below we list the methods and their properties. We use the notation 



to represent 
if a = 1, and 
to represent 



(aa) (22) 

(e aM e aM ...e aAN ) a (23) 
(aaf (24) 



(e aM e aM ...e aAN ) a (25) 
if a = — 1. So, for example, the 2nd order method 

(e M e M . ..e AN ) (e Ajv . . .e M e M ) = (e M e M . ..e A ») (f^e~ M . ..e~ AN \ l (26) 
is represented by 

(l)(lf (27) 

Note that the transpose of any method gives another equivalent method, as does 
permuting the entries in the fundamental unit. 

For odd order methods, the residue has an odd number of brackets in the commuta- 
tors. So, because the transpose of an individual bracket is minus that bracket, 

(odd order method) (same odd order method transpose) (28) 

gives a method of one order higher. For example, we can make a 4th order method from 
a 3rd order method, or a 6th order method from a 5th order method. 



5.1 Integer solutions 

The 3rd order integer methods that we have selected using the criteria of Section [| are 
given below. 





3rd Order Methods 


z 2 
z 3 

^3 

Z A 

^3 
^3 


(l) T (l)(l)(l)(l) r (-2)' r (l)(l)(l) 

(lf(4)(2)(-5f(2f(3)(2)(2f(l) 

(lf(2)(2)(-3f(lf(2)(lf 

(3)(-4) T (l)(3)(2f(l) 

(5f(7)(12)(-13f(l) 
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D 


L 


I 


L/D 


R/D 


Z 




6 


10 


9 


1.67 


0.2 


0.9 


zl 


12 


22 


9 


1.83 


0.6 


0.6 


zl 


6 


12 


7 


2.00 


0.4 


0.9 


Z 4 


6 


14 


6 


2.33 


1.7 


1.2 


Z 5 


12 


38 


5 


3.17 


98.8 


1.9 



And the 4th order integer methods are 





4th Order Methods 


Z 2 

z 3 
z 4 


(1)- J -(1)(1)' J, (-2)(1) T (1) T (1) T (1) T (1)(1) T (1)(1)(1)(1)(-2)' J, (1)(1) T (1) 
(lf(2)(lf(-3f(2)(2)(l)(2)W(-3)(2f(l)(l)(ir 
(lf(2)(3f(lf(-4)(3f(3)(-4f(l)(3)(2f(l) 
(6) T (-7)(lf(l)(5f(5)(lf(l)(-7f(6) 





D 


L 


I 


L/D 


R/D 


Z 


Z 1 


12 


20 


18 


1.67 


0.6 


1.3 


z 2 

^4 


12 


24 


14 


2.00 


0.8 


1.1 




12 


28 


12 


2.33 


4.6 


1.5 


z\ 


12 


40 


10 


3.33 


50.2 


2.2 



5.2 Irrational solutions 

The equations that we have derived can be solved for irrational solutions. We have been 
able to find the shortest 3rd order method analytically. It can be proven to be unique. 
It is 

(ai)(-a 2 ) T (-a 3 ) T (a 4 ) (29) 

where 

ai = 1 

a 2 = ~(b - Vl3 + 2 V5 + 2>/l3) 
a 3 = 1/ (1 + a 2 ) 

a 4 = -a 2 (1 + a 2 ) / (3 + 2a 2 ) (30) 

Renormalising to give a) = 1, we have method 

ai = 0.451525513208585723409578820 

a 2 = -0.630880954030002500791663663 ( 

a 3 = -1.136710925213995714728206549 ^ > 

a 4 = -1.219117392452583938929449032 
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accurate to 27 decimal places. This method has Z = 1.7. 

From this 3rd order method, we can generate the 4th order method 

(ai) (-a 2 ) T (-a 3 ) T (a 4 ) (a 4 ) T (-a 3 ) (-a 2 ) (ai) T 



(32) 



We have also found short fourth order methods. We assume a solution of symmetric 



form, using the ansatz aj-i+i 
the equations 



-a; and ai-i+i = —a*. For 1 = 6, this leaves us with 



i=l 
3 

i=l 



1 

2 



]T [a? + 20^0? - a\ 



(33) 
(34) 
(35) 



i=l 



to solve. 

Combining equations and setting ol\ 



1, we find solutions of the form 



a 2 

«3 



2 (a 2 x + a 3 y + 1) 



1/3 



(36) 
(37) 



where 

y = -a 3 (a 2 x 3 + 1' 

and x has four possible values depending on the a^s. From our ansatz, a 4 = — a 3 , 
«5 = — «2, «6 — — «i, «4 = — ds, a?, = —a 2 and a§ = —a\. 
For a 2 = —a 3 = —1, x = —1 giving the method 1Z\ 



a i 

«3 



l(2 + V2) 
-I (2 + 
-± (1 + v/2) 



0.675603595979828817023843904 
-0.675603595979828817023843904 
-0.851207191959657634047687809 



(38) 



This method has been found previously by Yoshida |]TTJ in the two-operator case, we see 
here that it is also a method for an arbitrary sum of non-commuting operators. This 
method has Z = 2.67. 

For a 2 = a 3 = —1, x is the solution of 

x 5 + 3x 4 + 3x 3 - 3x - 3 = 



giving the method 1Z\ 



0l = -1.075035037431900314780251056 
a 2 = -1.024607977441460486144230714 
a 3 = -0.550427059990439828636020342 



(39) 



(40) 
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This method has Z = 2.53, and is slightly better than the above method of Yoshida. 
For a 2 = —0:3 = 1, x is the solution of 



2a; 5 + 3a; 3 + 3a; 2 + 3 = (41) 

giving the method TZ\ 

ai = 0.938925888779098070854126976 

a 2 = -1.002122279211397565598116356 (42) 
a 3 = -0.563196390432299494743989380 

This method has Z = 3.56. 

And finally, for a 2 = a 3 = 1, x is the solution of 

x 9 + 3a; 7 + x 6 + 3a; 5 + 3x 4 + 3x 2 + 1 = (43) 

giving the method TZ\ 

ai = 1.087752928204421689142747144 

a 2 = -1.131212302433601022822197398 (44) 
a 3 = 0.543459374229179333679450254 

This method has Z = 4.39. 

We also searched numerically for other irrational solutions and found no short asym- 
metric solutions (i.e. shorter than the symmetric solutions found analytically). 



6 An Efficient Technique for Deriving Sub-optimal 
Higher Order Methods 

The technique for finding higher order methods described above used a first order method 
as a fundamental unit. We can also use higher order methods as fundamental units. This 
makes it easier to derive very high order methods, but the methods will be sub-optimal in 
the sense that we only generate a restricted set of solutions, which is unlikely to contain 
the method that is optimal with respect to any given criteria. 

The technique of using higher order fundamental units works as follows: 
The method of order o from which we form the fundamental unit is 



(45) 



and the fundamental unit is 



1 1 ^ / N \ 

n (e M . . . e^Y 1 = exp (3ba} £ A N + (3b° +1 r 

i=l J V n=l / 



(46) 
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Combining a succession j — 1, . . . , J of fundamental units with parameters bj and (3j 
gives 



j 

n 



i=l 



N 



exp EAV/E^v + E/¥: 



Therefore, to obtain a method of order o + 1, we require 



(47) 



and 



This technique can be iterated to get arbitrarily high order methods. 
As an example, we start with the first order method 

(1) 

by transposing, we get the second order method 



now, we solve the equations 



and 



E/?A>o 



E/¥? = ° 



(48) 



(49) 



(50) 



(51) 



(52) 



(53) 



A simple solution to Eqs. ([5^ and [5^) is 2 3 = l 3 x 8. Ordering is not dictated by 
the solution, so we choose the method which is its own transpose and hence 4th order 
accurate 

r m-i d r m-i r mi A 

(54) 



[(l)(l) r ]'[(-2)(-2) J j [(1)(1) J 
Again, we solve the equations 



K 



E ikCk > o 

k=l 



and 



A" 



E 7.4 = o 
fc=i 



(55) 



(56) 
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which have the simple solution 2 5 = l 5 x 32. Again, choosing the ordering so that the 
method is its own transpose, gives the 6th order method 



T 



where I = 594. 



'D(i) 



16 



(-2) (-2) T ' (4) (4) J (-2) (-2) 



(i)(iy ' (-2)(-2f (i)(i) 



16 



(57) 



7 4th and 5th Order Formulae for exp([A[,A2 



As a byproduct of our analysis, we can also use Eqs. (10|-12) to search for approximations 
to gates involving commutators. To do this, we set aj > and of = for X ^ 2. An 
approximation for a gate involving a commutator may be useful if only a subset of 
the generators of a particular group is available in hardware, but a given algorithm 
needs another generator of the group. For instance, if exp (— ia x At) and exp (—ia y At) 
are available in hardware, but exp (— ia z At) is not, then we need a way to generate 
exp (-| [a x ,(Ty] At). 

After some searching, we have been able to find one method for exp [A, B] to fourth 
order. It is 



2)' J '(2)-[(-l)(l)] 12 [(!)(-!)] 



(58) 



with residuals 



Pl2 


P11112 


P21112 


P11221 


P22112 


P12221 


P22221 


12.0 


1.0 


2.0 


0.0 


0.0 


-2.0 


-1.0 



This method can be combined with its transpose to give a 5th order method. 

8 A Simple Application 

To illustrate our methods, we have applied first, second, third and fourth order methods 
to the exactly soluble operator 



-iAt(o- x +a y +rT z ) 



COS 



(^3 At) - ^ sin (v^ At) 

-(i- 1) -75 sin (V3At) cos(V3At)+ii sin (\/3At 



(i + 1)^5 sin [y/3 At 



We used the first order method 



■iAtcry ^ — iAta z 



(59) 
(60) 
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1 2 

log(t) 



Figure 1: Here, we plot log(error) vs. log(time). Error is calculated according to Eq. (|64|). 
The lines from top to bottom correspond to the 1st, 2nd, 3rd and 4th order methods of 
Eqs. (H), ©, H) and ©. 



the second order method 

Z\ = (l)(lf = (e-^°* e -^°y e -^^ ^ e -iAta Ze -iAta Ve -iAta x ^ ^ 

the third order method 

Z\ = (l) r (l)(l)(l)(l) r (-2) r (l)(l)(l) 

g — iAta z ^ — iAta y ^—iAta x \ f^ — iAta x ^ — iAta y ^ — iAta z \ f^~iAta x ^ — %At<j y ^~iAta z ^ 

iAta x g~iAt<jy ^—iAta z \ f^—iAta z ^—iAta y ^—iAta x \ f^2iAtcr z ^2iAta y ^2iAta x ^ (62) 
g— iAta x g— iAta y ^— iAta z \ ( ^— iAta x ^— iAta y ^— iAta z \ f ^—%Ata x ^—%Ata y ^—iAta 2 

and similarly for the fourth order method 

Z\ = (l) T (l)(l) T (-2)(lf (If (If (If (l)(lf (l)(l)(l)(l)(-2f (l)(lf (1) (63) 

As a measure of the error, we take the difference between the a x , a y and o z compo- 
nents of the exact solution and our methods, Aa x , Aa y and Aa z . We then calculate the 
error 



E = ^(Aa x ) z + (Aa y ) z + (Aa z ) z (64) 

In Fig. ([]]), we plot the logarithm of the error as a function of the logarithm of the 
time that the system was evolved for. The first order method results are uppermost 
and higher order results lie underneath each other with fourth order results being the 
lowermost plotted. At = 0.01 for all methods. 
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Notice that the first order error oscillates once it reaches order 1. The rest of the 
errors remain small throughout the simulation, with the fourth order error remaining 
below 10~ 3 for the entire evolution. 

The error for all methods goes as nR(At)° +1 , where n is the number of times the 
method has been applied. Therefore, log-E = log R(At)° +1 + logn. For At = 0.01, 
this makes the y-intercept decrease by order —2 as the order of the method increases. 
Since the time evolved is proportional to n, the slope of the errors is 1 for all methods. 



9 Symplectic Methods 

In the study of classical Hamiltonian systems, we can cast the evolution of the coordinates 
qi and momenta Pi of fields or particles in the same language as we have done above for 
quantum systems. 

Write z = {qi,Pi)- Then the Hamilton equations for the system are 

z = {z,H} (65) 

where {a, b} is a Poisson bracket. Now, define D^z = {z,H}. The Hamilton equations 
become 

z = D H z (66) 
The formal solution to these equations is then, 

z(t) = e D ^z (67) 

Often, Dh can be separated into kinetic and potential parts Dh = Dr + Dy. In this 
case, we have the formal solution 

z (t) = e {DK+Dv)t z (68) 

Typically, symplectic methods approximate the above case (|B~8|) , in which there are only 
two operators in the exponential. Symplectic methods for two operators exist up to 8th 



order in the expansion ] I 



In our work, we have developed methods to approximate the case where there are an 
arbitrary number of operators in the exponential. This is important for simulations on 
both quantum and classical computers, since there can often be more than two terms 
which do not commute in the Hamiltonian. 

For example, any Hamiltonian of the form 

H = gij (q)PiPi + V(q) (69) 

where gij and V are functions of the q^s, can have an arbitrary number of terms which 
do not commute with each other. 
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A simple example of a quantum system where extra terms in the sum are necessary 
is an Ising spin system with next-nearest neighbor interactions. Here, the Hamiltonian 
becomes 

N 

H = 2 O 7 * ' CFi + 1 + a i ' °~*+ 2 ) ( 70 ) 
i=l 

In this Hamiltonian, none of the terms 0V2 • &i, 0V1 • cr?, 0$ • <7i+i or <7j • a i+2 commute. 
Therefore, for this system, we can arrange the Hamiltonian to have, at best, four terms 
which do not commute with each other. 



10 Conclusions 

The object of this paper has been to provide higher order approximation methods for op- 
erators of the form exp Y^Li An in terms of operators of the form exp (A±), exp (A 2 ), . . ., 
exp (An). We have focused on approximation methods of this kind since they are partic- 
ularly useful in quantum many-particle simulations for which the discretised Hamiltonian 
on a quantum computer takes the form of an exponential of a sum of non-commuting 
terms. 

To find higher order methods, we have derived and solved equations for methods up 
to 4th order. We find that the equations give a large number of methods, so we have 
selected a small number of them based on what seem to us to be reasonable criteria and 
presented them above. 

As a by-product of our search, we have also been able to find higher order approx- 
imation methods for operators of the form exp [A, B] in terms of operators of the form 
exp (A) and exp (B). These may be useful for quantum gates where exp (A) and exp (B) 
are available in hardware, but the gate exp [A, B] is desired for some particular algorithm. 

Our analysis has also shown that there is a quick technique for deriving approximation 
methods to arbitrarily high order involving the solution of relatively simple equations 
at each order. We have also presented these results, but it turns out that they lead to 
approximations that are far from optimal in the sense that there are many more gates in 
these methods than should be necessary. That is, they are accurate to high order, but 
relatively costly to implement. 

As an example of how useful our approximations can be, let us consider a case in 
which we want to apply an approximation method for time T = 1 with total error 
E = 1CT 4 . For a first order method, this means that we require about 5000 applications 
of the method. For second order, we require about 30 applications. For our third order 
method Z%, we need 2 applications. And for fourth order method Z\, one application of 
the method is more than sufficient. This results in a reduction of orders of magnitude in 
the computational cost of a given simulation or gate application. 

Using our equations, it is possible to search for 5th order methods (and from these, 
via transposition, to obtain 6th order methods). We made a number of attempts at 
the search, but were unable to find any 5th order methods due to the large size of the 
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search space. Thus, the only methods of 5th order and higher that we found were those 
methods mentioned above which tend to involve unnecessarily large numbers of gates. 
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Appendices 

A.l Proof of lower bounds on integral method sizes 



i=i 

7 

= E [{ a i-i + ^ a P ~ (^-0 



in 



i=i 



E 

i=i 



p-i 
E 



p! 



1 ?!(p-?)! 



)*(»!-,)"+« 



(71) 



ctj = ±1. Therefore, if p is odd, then 



E "fri = E = a i 
i=i i=i 



(72) 



and if p is even, then 



E aP i aP i = E aP = E (! - a i) aP + E = E ( x - "0 aP i + °j 



(73) 



i=l 



i=l 



1=1 



i=l 



Taking p = 2, the factor ^—yy , g = 1, is equal to 2, and the factor (1 — otj) is or 

2. A 2nd-order method requires aj = 0, therefore (a}) 2 must be even, and so aj must 
also be even. 

Taking p — 3, the factors ql ^i q y , q — 1,2, are equal to 3. A 3rd-order method requires 

<rf = 0, therefore (<r}) 3 must be a multiple of 3, and so aj must also be a multiple of 3. 
Taking p — 4, the factor ;^~yp 9 — 1) 2, 3, is even, and the factor (1 — a^) is or 2. 

A 4th-order method requires aj = 0, therefore (o"}) 4 must be even, and so aj must also 
be even. 

Taking p — 5, the factors q ^pl q y_ , 9 = 1,2, 3, 4, are multiples of 5. A 5th-order method 

requires erf = 0, therefore (aj) 5 must be a multiple of 5, and so aj must also be a multiple 
of 5. 

Combining these aj must be a multiple of 2 in a 2nd-order method, a multiple of 6 
in a 3rd or 4th order method, and a multiple of 30 in a 5th-order method. 
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A. 2 Formulae for the 

and A2 

The B 2 are defined by 



B^s in terms of the commutators 



e aAl e aM =ex P 5> p J Bf 
P =i 



The Campbell-Baker-Hausdorff formula, Eq. (f|), then gives 



Of Ai 



(74) 



B\ 2 



B 



13 



5 



14 



B 2 2 



B 2 , B 2 



R 1 R 12 



-^2) -^2 



R 2 R 3 
#2> -°2 



1 , 

- 2 A 12 

— (A U2 + A 221 ) 

— {A112 — A 22 \) 
1 

— 221 
24 

— (^1112 + A 222 i) 

2 (A1.112 — ^2221 — 2A1221) 
(A1.1112 — 2A21112 — 6A1122I 



720 



22112 



2A 



12221 



22221 , 



1 



(^4 



11221 



.4 



22112 1 



24 

2^(^12,^112] + [A 12 ,A 221 ]) 



(75) 
(76) 

(77) 

(78) 

(79) 

(80) 

(81) 

(82) 
(83) 



B^ ee 



B™ ee 



R 1 Rl3 
B 2 ,B 2 



r)2 r>12 
^2,^2 



(^21112 + A1122I — A 22 \\ 2 — Ai222l' 



(^4ini2 + A 2 \\\ 2 + A\ 222 \ + A 2222 i) 
- [A 12 , A ll2 } + [A 12 ,A 221 ]) 



12 

;< 

- (^21112 + ^4.11221 + A 22 \\ 2 + A1222I, 



(84) 
(85) 

(86) 



Rl d112 
-D2 5 -°2 



(^4 



11112 



+ A 



21112 



2Ai 



1221 



+ 2,4 



22112 



A 



12221 



-A 



22221 , 



(87) 
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A. 3 A simple measure of the error 

The error for a given method is given by Eq. ( |i~5|) 

r = £ of B* (88) 
x 

where X G {4, 13, 112} for a 3rd order method and X G {5, 14, 23, 113, 221, 1112} for a 
4th order method. 

r is a vector in the vector space of the commutators for which we do not know the 
metric. We would like to have a scalar measure of the error, and thus must pick some 
basis for the vector space. We choose the basis to be the commutators of A\ and A 2 . 
This basis is simple and spans the vector space of the -B/v's without redundancy. Since 
this basis spans the space of the B^'s, we do not need to go to iV larger than 2. For 
iV = 2, we can re-express r as 

t = Y,PyA y (89) 

Y 

where 

Y G {1112, 1221, 2221} (90) 

for a 3rd order method and 



Y G {11112, 21112, 11221, 22112, 12221, 22221} (91) 

for a 4th order method. The formulae for the py's in terms of the <xps are given in 
Appendix A. 4. In this basis, our measure of the error then becomes 

R= JEM 2 (92) 
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A. 4 Formulae for the py's in terms of the crps 

For N = 2, 



r = ]>>f = £p y A y (93) 

X Y 



Therefore, using the formulae in Appendix A. 2, we obtain 



Pi = °\ (94) 

92 = a) (95) 

P12 = \<J} (96) 

P112 = + (97) 

P221 = (98) 

Pni2 = ^ + ^ 12 (99) 

P1221 = y a ° a i-°Y 2 ( 10 °) 

P2221 = \-^f~\^ 2 ( 101 ) 
1 1 1 

1 „5 | 1 113 | 1 J112 /i no i 

P11112 = -^+12^ +2 a/ (1 ° 2) 

^5 ^ 23 | 1 113 | ^ 221 | ^ 1112 / 1 riQA 

P21112 = ^ - -a, + -a, + i<7/ + -a, (103) 

1 ^5 | 1 14 1 23 | ^ 221 ^1112 nn/A 

P11221 = + -a, - + -a, - a, (104) 

p 22 m = — o\ - -of + -a] 3 + -of 1 + a} 112 (105) 

H22L12 120 24 7 24 7 4 7 7 y J 

1 ^5 i ^ ^23 , ^ _113 , 1 221 ^ 1112 ( -\ r\a\ 

P1222! = ago*/ + 24*/ + 12*' 4^ ~ 2 7 (1 ° 6) 
1 1 1 

5 I ^.113 ^1112 

P2222! = -— <7 J + -<7, (107) 
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A. 5 Tables of Residual Errors 

3rd Order Integer Methods 





Pi 


P1112 


P1221 


P2221 


z 1 

^3 


6.0 


-1.0 


0.5 


0.0 




12.0 


-4.0 


-3.0 


5.0 




6.0 


-2.0 


1.5 


1.0 


21 


6.0 


0.0 


4.5 


9.0 


2\ 


12.0 


-864.0 


792.0 


180.0 





P11112 


P21112 


P11221 


P22112 


P12221 


P22221 


^3 


2.2 


3.1 


-3.2 


5.3 


0.1 


-1.3 


z 2 

^3 


13.4 


104.2 


105.6 


26.1 


84.2 


28.9 


z :i 

^3 


0.7 


5.1 


3.3 


1.8 


3.1 


1.2 


^3 


2.7 


8.1 


-2.7 


10.8 


-6.9 


-13.8 


■^3 


-3801.6 


-1900.8 


2505.6 


-1166.4 


499.2 


206.4 



3rd Order Irrational Method 





Pi 


P1112 


P1221 


P2221 




1.0 


0.012008 


-0.052816 


-0.058414 





P11112 


P21112 


P11221 


P22112 


P12221 


P22221 




0.001754 


0.003500 


-0.009304 


0.017412 


-0.014311 


-0.026310 



4th Order Integer Methods 





pi 


P11112 


P21112 


P11221 


P22112 


P12221 


P22221 




12.0 


-1.6 


0.2 


-3.4 


5.6 


-1.8 


-2.6 




12.0 


3.4 


6.2 


3.6 


3.6 


2.2 


-4.6 




12.0 


26.4 


40.2 


-5.4 


21.6 


16.2 


5.4 


n\ 


12.0 


-369.6 


-220.8 


309.6 


-86.4 


259.2 


86.4 



4th Order Irrational Methods 





pi 


P11112 


P21112 


P11221 


P22112 


P12221 


P22221 


n\ 


1.0 


-0.000414 


-0.008682 


-0.007027 


-0.026045 


-0.026732 


-0.004684 


n\ 


1.0 


-0.022171 


-0.013256 


0.014902 


-0.009176 


0.002796 


0.001717 


n 


1.0 


-0.001297 


0.038072 


0.035227 


-0.080082 


-0.079215 


0.001270 


n\ 


1.0 


0.002074 


0.196582 


0.194095 


-0.052861 


-0.050727 


-0.002155 
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